clear all
clc
close all

load basicgame

n=length(y);


%% Make parameter grid

tau=0:0.1:1.5;
lambda=0:0.01:1;

%% Run cross-validation

K=10;
c = cvpartition(n,'KFold',K);
levR{1}=1/3*ones(1,3);
levC{1}=1/3*ones(1,3);
besttau=cell(1,K);
bestlambda=cell(1,K);
ErrTL=zeros(1,K);
ErrLogitPCHM=zeros(1,K);
VarDiffLPCHM=zeros(1,K);
VarDiffU=zeros(1,K);
Covar=zeros(1,K);
empdistr=cell(1,K);

for k=1:K
    
    % Training Data
    
    train=find(c.training(k));
    trainy=y(train,:);
    idTrain=id(train);

    [besttau{k},bestlambda{k}]=CHlogit(tau,lambda,Row,Column,trainy,idTrain);
    
    % Table Lookup
    
    idTrain=id(train);
    r=length(unique(id));
    empdistr{k}=zeros(r,3);
    if r~=466
        'WARNING'
    else
        for i=1:r
            gamesemp=find(idTrain==i);
            actemp=trainy(gamesemp);
            empdistr{k}(i,:)=[length(find(actemp=='1')),length(find(actemp=='2')),length(find(actemp=='3'))];
            empdistr{k}(i,:)=empdistr{k}(i,:)/sum(empdistr{k}(i,:));
            clear gamesemp
        end
    end     
    
    % Test Data
    
    test=find(c.test(k));
    testy=y(test,:);
    idTest=id(test);
    
    T=besttau{k};
    L=bestlambda{k};
    [ErrTL(k),ErrLogitPCHM(k),VarDiffLPCHM(k),VarDiffU(k),Covar(k)]=combinederr(tau(T(1)),lambda(L(1)),Row,Column,testy,idTest,empdistr{k});
end

'Best Parameters'
besttau{1} %vary to see best-fit estimates of tau in different folds of CV
bestlambda{1}

'LogitPCHM Error'
mean(ErrLogitPCHM) 

%% Parameter Estimates

for k=1:K
    fitlambda=lambda(bestlambda{k});
    fittau=tau(besttau{k});
end

'Tau Estimate'
mean(fittau)

'Lambda Estimate'
mean(fitlambda)

'Table Lookup Error'
mean(ErrTL)

'Completeness'
kappa=(-log(1/3)-mean(ErrLogitPCHM))/(-log(1/3)-mean(ErrTL))

'Standard Error, Completeness'
se=sqrt((mean(VarDiffLPCHM) - 2*kappa*mean(Covar) + kappa^2*mean(VarDiffU)) / (-log(1/3)-mean(ErrTL))^2) * (1/sqrt(n))

save output fittau fitlambda ErrTL ErrLogitPCHM VarDiffLPCHM VarDiffU Covar  

    

